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NOTATION 


A parameter appearing in equation (B3) 

b impact parameter 

B parameter associated with detailed balancing (eqs. (B5) and (B6)) 
c. molecular speed 

1 /2kT" 

c G most probable molecule speed, J/' — — — 

c 5 v most probable molecule speed at the initial temperature T]_. 

C exponent in selection rule (eq. (7)) 

d strength associated with inter molecular potential 
E kinetic energy 

E r rotational energy, k0 r j ( j + 1) 
f distribution function 

g relative speed 

I moment of inertia 

j rotational energy level 

j D average energy level where 0 r j o (j o + 1) =* T e 
k Boltzmann constant 

m magnetic quantum number; also molecular weight 

Mj. momentum, /2lE r 

N number of simulated molecules within the cell 

n number density 

P rotational transition probability 

P modified rotational transition probability 

Q collision cross section 

R intermolecular distance 
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t time 

T temperature 

u x component of velocity 

v y component of velocity 

V interaction potential, or velocity 

w z component of velocity 

6 index of power associated with point centers of repulsion model, 
(eq. (6)), V (R) = d/R<$. 

A rotational level jump 

© r characteristic rotational temperature 

A relaxation time, (eq. (16)) 

p reduced mass 

v collision frequency 

a effective collision diameter 

x characteristic collision time based on equilibrium translational 
temperature, l/(n-rra^c 0 ) 

X deflection angle defined in equations (2) or (10) . 

Q solid angle 

Subscripts: 
e equilibrium 

i,j rotational states 

max maximum value 

r rotational 

t translational 
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Superscript 
vector 
r after 

U) £ th momentum 
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SUMMARY 


Theoretical studies of translational and rotational energy relaxation in 
diatomic gases are described. The direct simulation Monte Carlo method is 
employed to solve the Boltzmann equation for a rotationally excited highly 
nonequilibrium gas. The gas investigated is homonuclear diatomic nitrogen, 
and the semiclassical model of Itikawa is incorporated for the transition 
probability that describes rotation-translation energy interchange. 

The details of energy interchange between the translational motion and 
the rotational energy levels of the gas are examined for spatially uniform 
flow without boundary interactions (the "box" calculation) with a variety of 
initial conditions. The results show: 

1. The assumption that relaxation .occurs via- successive, local Maxwellian 
velocity distributions, which is a commonly used basis for finding approximate 
solutions of Boltzmann equation, is not valid for gases that are initially in 
highly nonequilibrium states. This is especially true for initial conditions 
that involve low translational and high rotational temperatures. 

2. The energy distributions for such transitions show bimodal (or double 
peak) relaxation patterns; the secondary peak ("satellite peak") appears around 
the Maxwellian elastic peak in the velocity distribution early during the 
relaxation period. The secondary peak is due to inelastic collisions and is 
analogous to the rotational Raman effect accompanying Rayleigh scattering. 

3. The rotational energy distribution also shows bimodal relaxation 

effects: In addition to thermal equilibrium Boltzmann peak, a weak peak also 

appears at the high rotational energy levels. When the rotational energy 
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distribution is a delta function, however, relaxation proceeds only as a 
single-peak distribution. One, therefore, concludes that single- or double- 
peak relaxation depends on the type of initial distributions assumed. 

4. Relaxation of the velocity distribution to equilibrium Maxwellian 
occurs relatively fast while the rotational energy relaxes more slowly. The 
relaxation time depends not only on equilibrium temperature, but also on 
initial velocity and rotational energy distributions. 

Close correlation of the relaxation between the box models and fluid 
flows, such as, sound absorption, shock wave, and free-jet expansion experi- 
ment are described. Also presented are brief preliminary results of a shock 
'wave showing translational and rotational energy relaxation structure. 

A 16-mm movie film displays examples of the relaxation effects of the 
"box" model with a variety of initially specified velocity and rotational 
energy distributions. 


INTRODUCTION 


[ 

A knowledge of internal energy transfer mechanisms at the molecular level 
is valuable for an accurate understanding of many important 1 nonequilibrium 
problems that occur in high-speed gas dynamics, acoustics, laser transmission, 
detonation, combustion, pollution, and atmospheric physics. For example, 
collision- induced rotational transitions play a major role in establishing the 
population inversions leading to gas-dynamic laser action, an d als o in evalu- 
ating the effects of highly nonequilibrium energy transfer in rarefied gas 
flow about spacecraft entering the planetary atmosphere. 

In the present paper, several new and important results are presented on 
internally excited translation-rotation energy relaxation. These results are 
obtained by solving the Boltzmann equation by the Monte Carlo direct simula- 
tion method, which previously has been applied successfully to monatomic 
(without internal rotational relaxation) gas flow problems (refs. 1 to 3). 

An important feature of the simulation method is that it provides insight into 
the effects of collisional relaxation at the microscopic level. In particular, 
the instantaneous internal energy distributions can be continuously observed 
throughout the relaxation processes. To ensure that these distributions be 
meaningful, however, it is essential that the rotational transition probability 
function, used in the method, display certain features; namely, (1) probability 
must be conserved, (2) probabilities relating transitions to and from pairs of 
definite states must satisfy "detailed balancing," and (3) probabilities, when 
used in the simulations, must yield the correct asymptotic behavior of the dis- 
tributions (refs. 4 and 5). 

The Monte Carlo method itself can be described briefly as follows (see 
refs. 2 and 3 for details): The flow is determined by following a statistical 

sample comprising several thousand molecules that are allowed to collide with 
each other. The phase space coordinates that involve trajectory and rotational 
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variables are known at every instant. These coordinates are allowed to change 
only during a collision; the modeling of these intermolecular encounters is, 
of course, the essential part of an accurate simulation. To account for these 
sampling encounters, a molecule and a near neighbor are each selected at 
random as are also their impact parameter and relative orientation angles — 
all in a manner representative of typical molecules undergoing similar 
encounters. They are accepted for a collision or rejected according to a 
selection rule that is dependent on the collision cross section. Since the 
initial coordinates, that is, relative velocity, impact parameter, orientation 
angles, and pair of rotational quantum states, are known, the final rotational 
quantum states can be computed. This involves computing the transition prob- 
abilities of all quantum states that are accessible from the known initial 
states and then selecting randomly from this resulting distribution. 

The procedure used here for the "translational" interactions parallels 
other investigations (refs. 1 and 3) which. treat of monatomic gases only. The 
procedure is different, however, from those investigations that have treated of 
translation-rotation interactions. All investigations are easily categorized 
under the following descriptions: (1) semiempirical, (2) classical, (3) semi- 

classical, and (4) quantum mechanical. Within these categories, the semi- 
empirical treatment includes an energy sink (ref. 6) and rough spheres and 
loaded spheres (ref. 7) to model the translation-rotation collision processes. 
While such methods do not appear to be satisfactory for highly nonequilibrium 
flow, they adequately describe near equilibrium steady flows. The classical 
models (refs. 8 and 9), although consistent with the classical direct simula- 
tion Monte Carlo procedure used here, necessarily include approximations to 
make the models tractable for studies of the type considered in this report. 

The approximations yield appropriate macroscopic behavior for a nonequilibrium 
example, but do not adequately provide limiting microscopic beh avior. In 
particular, individual molecular encounters that violate energy and momentum 
conservation can occur. 

Semiclassical methods (refs. 4, 5, and 10) appear to have physically 
realistic bases. The simplified model of Pearson and Hansen satisfies limit- 
ing equilibrium behavior, but, during a calculation, the model causes a drift 
in the answers that violates energy equipartition (ref. 4). Itikawa's model 
is more rigorously founded, allows for treatment of molecular collisions (ref. 
5), and also satisfies conservation of probability and appropriate detailed 
balancing. Itikawa's model, therefore, satisfies the desirable characteristics 
of the ideal model that we described earlier; our investigations described in 
this paper depend on this model. Our intent, then, is to extend its applica- 
tion to even more general problems. 

As regards the fourth category of the model (i.e., quantum mechanical 
models), the author is not aware that truly quantum mechanical results are yet 
viable. Such descriptions are difficult to obtain analytically, and to apply. 

In this paper, we treat translation-rotation interactions for a spatially 
uniform gas far removed from solid boundaries. We are concerned only with a 
basic understanding of translation-rotation relaxation behavior in highly 
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nonequilibrium situations. In’ fact, it is our belief that the Monte Carlo 
method is best suited for studies of the type described in this paper. 

The results presented are based on calculations involving three differ- 
ent basic types of initial conditions: (1) equilibrium, (2) nonequilibrium- 

equipartition (i.e., equipartition is satisfied, but distributions' are 
perturbed), and (3) nonequilibrium-nonequipartition (i.e., both the equiparti- 
tion and the distributions are perturbed). Also included are the results of 
monatomic gas simulations (rotational relaxation effects frozen) to assist com- 
parisons with coupled translation-rotation relaxation simulations.. To further 
assist the understanding of the Monte Carlo method, the essential mathematical 
relations are also given in this report. 


FORMULATION AND PROCEDURE 

The essence of the Monte Carlo procedure is described briefly in the 
Introduction. Introduced in this section are several analytical relations 
that assist both the understanding and use of the method. Appendices A and 
B provide supplementary analysis to the procedure and appendix C is a listing 
of the computer program that was used in the procedure. 


Governing Equations 

The study described in this paper concerns the temporal and spatial relaxa- 
tion of the velocity and rotational energy state distributions that character- 
ize a statistical sample representing several thousand molecules. If we assume 
that the molecular distributions themselves are diagonal and independent of the 
degenerate rotational m substates, the Boltzmann equation (or Wang, Chang 
and Uhlenbeck equation, see ref. 11) that relates the temporal and spatial 
behavior of the distribution functions f^ can be written 


9 
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where n is the number density, = f^ (x, V, t) is the distribution func- 
tion which depends on time t, position velocity and rotational state 
i, dcr/dft is the differential cross section corresponding to solid angle 
External forces are assumed to be absent. The Monte Carlo procedure is used 
to effectively solve this equation by means of a probabilistic sampling pro- 
cedure. Implicit within the equation and procedure are the conventional fluid 
dynamic conservation laws (i.e., conservation of mass, momentum, and energy; 
see, e.g., ref. 11). 


Of greatest interest for the study given here is the "box" calculation 
wherein the gas is spatially uniform, has constant density, and is stationary; 
that is, the gas is entirely contained within an "imaginary box" that has unit 
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volume and noninteracting boundaries (sketch (a)). The procedure, in this case, 
concerns interactions in a closed system. Energy conservation is applied 
directly (i.e. , exact energy conservation is imposed on the interacting pairs 
of molecules) and "random sampling" ensures that, over long periods of time, 
the number of molecules contained within the box remain constant. These con- 
cepts are treated in greater, detail in the subsequent discussion. 



The "box" calculation has general utility since such a caxcuranon, wnen 
started with an appropriate set of initial conditions, provides insight into 
mechanisms in the more general flow situation as found in sound absorption, 
normal shock wave, or free jet expansion experiments (see table I). In addi- 
tion to the "box" calculation results, a simulation is also given for steady 
one-dimensional shock wave flow. This result is preliminary and demonstrates 
the capability of the method for simulating more complicated flows. 


Collision Parameters 

The essence of an accurate simulation is the random or probabilistic 
sampling used to select the interacting (colliding) molecular pairs, to deter- 
mine whether a reaction occurs, to find the resulting "states," and then to 
advance the time interval for the next collision, and so on. To provide 
insight on this entire collision process and to arrive at a criterion for 
evaluating certain of the parameters required to define a collision, it is 
worthwhile to briefly review the classical representation of the equivalent 
process, and to observe how such relations depend on intermolecular potential. 

A classical representation (ref. 7, ch. 8) is given by 

CO 

X (b , g) = TT - 2 / (b dR/R 2 )/Vl - (b/R) 2 - V(R) / (1/2) yg 2 (2) 

R C 7 

r 00 

Q^(g) “ 2ir J (1 - cos^ x)b db (3) 

o 

where V(R) is a spherically symmetric intermolecular potential, X (b ,g) is 
the encounter deflection angle, which depends on impact parameter b and on 
g s relative velocity of approach, p is the reduced mass, R<. is the distance 
of closest approach, and Q^^(g) is the &th "momentum" transport cross sec- 
tion, which also depends on relative velocity (for studies reported here, 

Z = 1) , The collision frequency, v, is then given by 
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(4) 


v 



With this relation, the collision time, At, between encounters and the elapsed 
time, t, are given, respectively, by 


A 2 1 

At = — — 
N v 


£ - ? At i 
x 


(5) 


where N is the number of particles in the "box." 

Collision . — For cases in which the intermolecular potential is inversely pro- 
portional to the power, 6, of distance between colliding molecular pairs, we 
can readily calculate a frequency ratio (e.g., ref. 5) 



( 6 ) 


where V max and g max are the maximum possible values in a cell. The dimen- 
sionless ratios of frequency and relative approach velocity are related through 
6. If the collision process can be represented by symmetric inverse power-law 
potentials, then equation (6) is a valid representation for all collisions 
and can be used as a criterion to decide whether an "encounter" is a "collision. " 
The representation, therefore, serves as a "selection rule" for encounters. 1 

For the Monte Carlo results displayed in this report, we have arbitrarily 
picked intermolecular potentials with 5=4 (i.e., "Maxwell molecules"). 

Of all encounters that are collisions, we must further categorize those 
which are elastic from those which are inelastic (i.e., those which yield 
rotational transitions). 

Inelastic collision . — Not all collisions result in a rotational transition. 

For example, some interacting pairs may have insufficient relative energy to 


1 Actual intermolecular potentials have a more complex behavior than the 

idealized 6 potential upon which equation (6) is based. Relations equiv- 
alent to equation (6), but based on more accurate representations for the 
molecular potential, can be found in reference 12. The expressions were 
derived recently and, hence, were not available for the simulations described 
in this paper. 
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induce a transition. To separate such events from those which result in 
rotational transitions (i.e., inelastic collisions), we introduce a relation 
similar to that given above but with a different value for the exponent, 
that is. 


V 

max 1 . 

A 

where the subscript i denotes inelastic collision. The appropriate, value to 
use for C^, however, is based op simulation results. We require that simula- 
tions, which start with Maxwell-Boltzmann distributions that .satisfy equiparti- 
tion, must yield nondrifting results. A value of Ci = 0.431 yields this 
desired behavior. The procedure used to evaluate C-j^ is also described in 
reference 5, but in greater detail. 




Collision Dynamics 

In the previous section, the parameters required to determine the occur- 
rence of a collision are given. In this section, we describe 'the procedure 
for finding the trajectories after a collision. A collision, of course, also 
can be accompanied by rotational transitions in either or both colliding pairs 
of molecules. These rotational transitions can also perturb the particle 
trajectories. In this section, we describe the relations that ensure colli- 
sion symmetry (i.e., a collision is invariant with its inverse) and that enable 
rotational transitions to be more precisely determined. . 

The relative velocity and impact parameter after a collision are obtained 
by knowing the onset energy and momentum. The relations are given by 


and 


(«’)2 - g2 - - CE^ - E rl + e; 2 - E r? ) 


b' = 


[gb - (m; x - M rl + m; 2 - M r2 )/y] 


( 8 ) 

(9) 


where E r and M r are the rotational energy and momentum before a collision 
and a prime distinguishes the corresponding values after a collision. We 
assume that rotational transitions only slightly perturb the relative velocity 
and impact parameter. 

The deflection angle, given by equation (2) , is actually dependent on the 
functional behavior of the intermolecular potential, but, for finding the 
limiting trajectories, we assume that the infinite-rise potential (i.e., a 
"billiard ball" collision) is adequate. . The corresponding deflection angle is 
given by 
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X(b) = 2 cos 1 


( 10 ) 


where O is the effective diameter of the rigid-spherical molecule. The 
slightly perturbing effect of inelastic collisions resulting from the rota- 
tional transitions is accounted for by the following average: 


x<b) = x(b) + 2 x(b,) 


(id 


To completely specify a collision, however, it is also necessary to give the 
orientation angle e, which references the collision plane with respect to 
some arbitrary coordinate plane (see, e.g., ref. 13, p. 36). The velocity 
components before and after the collision can be related (e.g., ref. 13) 


g' = (g cos X - /g z - g z cos e sin 

X § X X 

s y " [ g y cos * + (s x g y cos e + 8 g z 
g 2 = \ [s z cos X + (g x g z cos e - g g y 


X) 


sin e) sin x/v'g* - g z 
sin s) sin x//g z - g 2 

X 


] 

] 


(12a) 

(12b) 


.(12c) 


We impose conservation of linear momentum to find the resulting velocity com- 
ponents after a collision. There results 



1 

2 

1 

2 

1 

2 


(Ui + u 2 + 

( V 1 + v 2 + 
(w^-+ W2 + 





(13a) 


(13b) 
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Rotational Transition Probability 


To describe inelastic collisions, also needed, in addition to the trajec- 
tory parameters introduced in the previous sections, are expressions which 
relate the probability of transition between initial and final rotational 
energy states (i.e., the rotational transition probabilities). As was pointed 
out in the introduction, the semiclassical probabilities derived by Itikawa 
(ref. 5) are used. A brief description of their properties follows. 

We describe an interaction where rotational transitions occur from levels 
(i,j) to (i ' , j ' ) by 


N 2 (i) +N 2 (j) •*N 2 (i’) +N 2 (j') 


(14) 


The collision trajectory itself, as described earlier, is determined classi- 
cally: Given an analytical relation for the trajectory, the amplitude of the 

rotational transitions can then be determined from quantum mechanical con- 
siderations. By appropriately combining the trajectory with an expanded set 
of Schroedinger equations (e.g., see ref. 5), the amplitude of the rotational 
transitions can be obtained by solving a set of coupled differential equations. 
In order to reduce the rank of the system, the effective potential method of 
Rabitz (ref. 14) is employed. The method eliminates the dependence of the 
interaction matrix on the magnetic quantum number, m. The resulting coupled 
set of ordinary differential equations are then solved by using the exponen- 
tial approximation (see ref, 15). What is important is that the method treats 
an interaction regardless of its "strength," and, in addition, allows for the 
likelihood of all transitions, including those with multilevel jumps. The 
simultaneous transitions for both colliding molecules (i.e., rotation-rotation 
as well as rotation-translation) are also taken into account. The precise 
formulation used is given in reference 5. 

Some important properties of the probabilities that pertain to the Monte 
Carlo simulation method are described briefly in appendix B. If a collision 
is inelastic as selected by equation (7), the transition of pair's molecular 
states are then determined by the Itikawa 's rotational transition probabilities. 


RESULTS AM) DISCUSSION 


In this section, Monte Carlo simulations are described for a stationary 
homogeneous molecular gas (i.e., for a "box" calculation). The simulations 
differ depending on the choice of the initial distribution functions (see 
table 1). The initial conditions fall under three general categories: • 

(1) complete equilibrium, (2) nonequilibrium, equipartition, and (3) non- 
equilibrium and nonequipartition. We use the term nonequilibrium here to 
denote that either the velocity distribution, ft, or the rotational energy 
distribution, f r , are non-Maxwellian. 
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The first case, that is, complete equilibrium, tests the method, as well 
as the internal parameters, for self-consistency. The velocity and r.o.tational 
energy distributions -should remain constant for extended periods of calcula- 
tion. That is, the internal energy distributions should remain Maxwellian and 
equipartition should be inviolate (e.g., refs. 4 and 5). 

The second case, where the velocity and rotational energy distributions 
are specified to satisfy energy partition (i.e., the fraction of energy dis- 
tributed between translation and rotation is proper, but where the distribu- 
tions, themselves, are non-Maxwell ian) provides a test on whether the procedure 
has an internal driving mechanism that will yield a relaxation to equilibrium 
within a physically realistic time. 

The third set of initial conditions, where the initial distributions 
violate both equipartition and are highly nonequilibrium states, allows -even 
more complex investigations. For example, one can study the relaxation 
processes to equipartition as well as how the velocity and rotational distri- 
butions interact during the relaxation. In effect, these simulations with 
varied starting conditions give qualitative information on the coupling of the 
energy distributions and quantitative data on the rates of relaxation. 

In table 1, the specific choice of initial conditions for the simulations 
described in this paper are listed. These results are also useful in providing 
qualitative information that can be used to interpret results in several 
equivalent experiments. The experiments are sound absorption, shock-waves, 
and free-jet expansions (ref. 11). Table 1 lists the simulation and the 
related experiment type. The simulations are described in the discussion that 
follows . 

Initial distributions; complete equilibrium. — The first test of a good method 
for simulating solutions to Boltzmann’s equation is that Maxwellian energy 
distributions, both in velocity and in rotational energy, not change for 
extended calculation periods. In figure la are given the results of such a 
simulation. The results show sets of paired figures for progressively increas- 
ing times corresponding to t/f = 0.0, 1.0, 5.0, and 10.0. One figure in the 
pair is a snapshot of the rotational energy distribution function, - f r , plotted 
versus rotational energy level j at a definite time t/x, and the other 
figure is the velocity distribution function, f t , plotted versus velo city 
c/cg, where eg is the most probable molecular speed defined by eg = /2 kT e /m, 
m is the molecular mass, and T e is the equilibrium temperature associated with 
the "box" model. We observe that, although small fluctuations occur around 
the dotted curves (which represent the true Maxwellian distributions) during 
the calculation period, these fluctuations do not grow (ref. 5). In fact, 
figure 2b shows the results of the same calculation, but where the time average 
of the distributions, given by 

I t -A/\ 


or t 


dt 


(15) 


10 



are plotted. We observe that the fluctuations are negligible in the second 
group of "snapshots." These simulations illustrate that, indeed, the pro- 
cedure is stable over long calculation periods. The next case tests the 
capability, of the procedure to drive arbitrarily specified initial distri- 
butions to the Maxwellian limit. 

Initial distributions: non-Maxwell ian velocity and rotational effects frozen . — 

The Monte Carlo method allows for considerable flexibility regarding the 
precise specification of the initial distributions. For example,- one can 
freeze the rotational relaxation effects and investigate only the relaxation 
of the velocity distribution. The next example is of this type. * In figure 2 
are displayed the resulting time history for relaxation of the velocity 
distribution function, starting with two different initial distributions. In 
figure 2a are displayed the relaxation process es th at correspond to initially 
letting every molecule have a speed equal to v^3/2 cq. The dotted curve is a 
Maxwellian distribution characterized by the temperature T e = 320 K. In 
this example, the "Dirac delta function" type of initial distribution should 
relax to coincidence with the dotted curve. The rotational energy, of course, 
is ignored. We observe that the distributions are largely equilibrated by 
the instant t/x = 1 (the area differences, between the solid and dotted 
c urve s correspond to the number of molecules that still have initial velocity 
1 / 3/2 cq and remain to be "equilibrated" — i.e., about 10 percent of the 
total). At t/x = 2.0, the distribution is established and very little change 
occurs thereafter. One concludes from this simulation that the procedure leads 
to the correct Maxwellian limit, as indeed it should. In figure 2b, the ini- 
tial distribution is slightly different. In this case, the energy, corre- 
sponding to kT e , is distributed at two separate initial velocities: co/2 and 

/Il/2 cq. The rotational energy is managed in the same manner as the example 
in figure 2a. The result for this case is nearly the same. In fact, little 
difference can he observed in a comparison of the relaxation history. At the 
instant t/r = 1.0, roughly the same fraction of molecules remain to be 
equilibrated as in the first example. The distribution appears to be estab- 
lished by the instant t/T = 2.0 and changes very little thereafter. 

Initial distributions: Maxwellian velocity and equipartition . — Our next 

simulation, figure 3, illustrates the relaxation effects that occur when the 
initial velocity distribution is Maxwellian and the rotational distribution, 
which satisfies energy equipartition, approximates a delta function centered 
at jo = 10 (i.e., every molecule has kT e rotational energy in the 10th 
energy level) . This rotational level also represents the probable rotational 
energy level for a Boltzmann distributions at temperature kT e (i.e., jo is 
found from k0 r jo(lO + 1) - kT e , where 0 r (N)2 — 2.9 K and T e = 320 K) . 

Since our investigation concerns homonuclear nitrogen, only rotational 
transitions that satisfy the multiples of Aj = ±2 are allowed. At the first 
instant displayed in figure 3 after relaxation begins (i.e., at t/x = 0.5), 
we observe a double peak appearing in the velocity distribution. This behavior 
is very similar to the Stokes and anti-Stokes lines that appear in Raman 
scattering (ref. 16). The position of these peaks can be calculated (see 
appendix A) and appear at the velocities c/cq = 0.91 and 1.1, respectively. 
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These peaks appear because the rotational energy is "dumped" into a narrow 
energy band. The transitions, Aj = ±2, are most probable compared to multi- 
level transitions, |Aj-j_>4. This effect results in two perturbations appearing 
in the velocity distribution about the most probable velocity cq. 

Initial distributions: Maxwellian velocity (Tt > Te) and nonequipartition . — 

The simulation displayed in figure 4 demonstrates the relaxation effects caused 
by violating equipartition. Rotational energy states are assumed initially to 
be unexcited; the energy is contained entirely within the velocity distribution 
which is Maxwellian (T t = 534 K) . The fact that single level (Aj = ±2 for 
homonuclear molecules) rotational transitions are the most likely compared to 
multilevel (|Aj[>4) transitions is also apparent here." After relaxation begins, 
we observe that the lowest level rotational states populate first. The gain in 
rotational energy appears to be at the expense of the molecules with velocities 
that correspond to the most probable velocity eg or higher. As the time t/x 
increases, energy continues to be "dumped" from the translational to the rota- 
tional mode as demonstrated by the downward drifting velocity distribution and 
the^ upward drifting rotational distribution’. These processes become less 
efficient as energy is "dumped" to higher and higher rotational energy levels. 
This is apparent because, as the width of the energy level increases, the 
energy interchange between the rotational and translational mode becomes even 
less efficient. We find, then, that considerable time is required to populate 
the uppermost rotational energy levels. In fact, the velocity distribution 
appears to be nearly equilibrated by the time t/x = 20.0, while the rotational 
distribution is still relaxing at t/x = 50.0. 

The simulation demonstrates that the step-wise populating mechanism 
implicit within the Itikawa model leads to relatively slow relaxation to a 
Boltzmann rotational distribution. Of course, if multiple level transitions 
were more effective, the rate of relaxation to a Boltzmann rotational distri- 
bution would be greatly enhanced. These features are characteristic of 
translation-rotation transitions, and they are apparent in all simulations 
involving translation-rotation interactions. 

The simulation displayed in the next figure, figure 5, differs from this 
example in that rather than follow the populating of the rotational energy 
levels from an initial state of "excessive" translational energy, we follow 
the depopulating of the rotational energy levels from an initial state of 
"excessive" rotational energy. Of course, "excessive" refers to the manner 
in which that energy is initially distributed relative to the distribution 
that satisfies equipartition. 

Initial distributions: Maxwellian velocity (Tt < T e ) and nonequipartition. — 

In figure 5a, we observe that, at the initial instant prior to relaxation, 
the rotational energy is stored in a Boltzmann distribution (T r = 793 K) which 
peaks near the most probable levels j = 10 or 12 (i.e. , f r (x) is maximum at 
x = 11.16). The velocity distribution, however, peaks at low velocities, 
c*/c 0 = 0.135. 

Several interesting features can be observed during the relaxation pro- 
cesses. We note that a satellite peak (refs. 16 and 17) develops on the high 
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sides of the velocity distribution (i.e., c/cg ~ 0.44), and this peak forms at 
the expense of rotational energy near the most probable levels j = 10 or 12 
(as exemplified by the "dip" in the rotational distribution) . This peak con- 
tinues to grow as rotational energy is converted to thermal motion. Again, one 
observes the effect of the inefficient coupling of translational and rotational 
energy interchange with the higher rotational energy levels. This effect, of 
course, results in the appearance of a second peak at the higher rotational 
levels. This effect of double peaks is also discussed by Polanyi and Woodall 
(ref. 18). The peaks remain in the distribution until quite late during the 
relaxation process (say, t/t = 50.0). We also notice that, by this time, the 
velocity distribution has nearly equilibrated to Maxwellian, and this occurs 
before the rotational distribution becomes Boltzmann, similar to what occurred 
in the previous example, figure 4. 

Also interesting is a comparison of how close the distributions approxi- 
mate the Maxwell-Boltzmann distribution during each instant of the relaxation 
processes. Such comparisons are displayed in figure 5b. Here, the dotted 
curves are "local Maxwell-Boltzmann" distributions rather than the asymptotic 
limiting equilibrium distributions displayed heretofore. The dotted distri- 
butions are determined by matching the energy, in both the velocity and 
rotational modes, with the simulation results. These results illustrate that 
the actual distributions found by the simulation deviate significantly from 
local Maxwell-Boltzmann distributions. This demonstrates that the popular 
methods that rely on expansion procedures involving local Boltzmann distribu- 
tions for solving the Boltzmann equation can be unreliable. In fact, the 
double-peaked results displayed in figure 5 illustrate that appropriate distri- 
butions can have rather complex non-Boltzmann functional behavior. . 

Initial distributions: Maxwellian velocity (Tt < T e ) and nonequipartition . — 

The simulation displayed in figure 6 is very similar to that displayed in 
figure 3; the initial rotational energy distribution approximates, a Dirac 
delta-function, but the simulation differs in that the constant rotational 
energy assigned each molecule violates equipartition. Here we have an initial 
dumping of the rotational energy into the 16th rotational energy level. The 
velocity distribution, however; is Maxwellian. At the onset of relaxation, 
we observe the "anti-Stokes" Raman scattering effect appearing in the velocity 
distribution, that is, the appearance of a satellite peak (using the relations 
in appendix A with j = 16 and Aj = -2, we find the' peak location to be 
c/cg = 0.55). The pumping mechanism by which the high rotational energy is 
converted into thermal motion also appears conspicuously in this example (see 
ref. 17). The single step transitions from the 16th to 14th level (Aj = -2) 
occur first with the quanta of rotational energy being preferentially absorbed 
by the very low-speed molecules. As this process proceeds, the number of 
molecules in the 14th level approaches equality with the number in the 16th 
level. Simultaneously, the "anti-Stokes Raman" peak broadens as the slow 
speed molecules are "pumped" into this region. As the relaxation progresses, 
the lower lying rotational levels are successively populated while the "anti- 
Stolces Raman" peak becomes increasingly broad. 

Contrary to the previous example, figure 5a, the rotational distribution 
relaxes continuously with one one peak. The velocity distribution, however, 
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similar to that in figure 5, displays a slow relaxation to the Maxwellian 
distribution. The relatively large energy in the high rotational states, the 
slow relaxation experienced in these states., and the strong coupling of these 
levels with the slow molecules all appear to contribute to inhabiting the 
relaxation of the velocity distribution to Maxwellian. 

To show the manner by which the local distributions compare with local 
Maxwell-Boltzmann distributions during the relaxation processes, the dotted 
curves are again introduced in figure 6b, as was done in figure 5b. Again, 
we observe that the distributions do not approximate Maxwell-Boltzmann distri- 
butions during each instant of the relaxation. 


Relaxation of Average Rotational Energy 

- The results displayed in figures 1 through 6 illustrate, in particular, 
the manner by which an initial energy distribution relaxes to the final 
Maxwellian distributions. Also interesting is the manner by which the "average 
energy" approaches some asymptotic constant value. Such results are given in 
figure 7. In this figure are displayed four curves. Three curves display the 
average energy relaxation associated with initial distributions, which are 
Dirac delta-function type, and the fourth displays results with an initial 
distribution that is Boltzmann with high rotational temperature. The essential 
feature is the comparison of results between curves that have high and low 
initial rotational energy (e.g., cases 4 and 5). Because the coupling between 
translation and rotation is efficient for the lower levels, the slope of these 
curves is greatest. Also, interesting is case 3 which illustrates that even 
though the initial distribution satisfies energy equipartition, the system is 
not bound to satisfy equipartition during the subsequent instants as the rota- 
tional distribution asymptotically approaches a Boltzmann distribution. Addi- 
tional simulations all illustrate the downward shift (as illustrated in the 
figure for ‘j = 10) followed by the upward relaxation to "equipartition." 

Figure 7 can also be used to define a useful relaxation time that charac- 
terizes the simulation results.- Such a definition, of course, is not exactly 
clear because the relaxing curves are not exponential. One can resort to the 
definition given by (e.g., ref. 11) 

A = [ (E r ) e - E r (t)]/(dE r /dt) (16) 


This definition, however, is impractical when the energy difference in the 
bracket is small. One can also define the relaxation time to be that time when 
the bracket expression has reduced to 1/e of its initial value (e.g., ref. 13). 
On the basis of this latter definition, it turns out that A — 32 is satis- 
factory for both curves labeled j = 0 and j = 16. This value also seems to 
be consistent with the simulation displayed for all three initial delta func- 
tions of rotational energy distribution (i.e. , j = 0, 10, and 16) in figures 
3, 4, and 6. In these figures, we can see that the relaxation appears to be 
nearly ceased at the same .instant between the two displays of the distribution 
function at t/x = 20 and 40. 
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The relaxation of translational energy is not shown here since, for our 
"box" calculation, conservation of energy E t is directly related to the 
average rotational energy E r via the relation 

E t /(E t ) e = 5/3 - 2/3 E r /(E r ) e (17) 


at each instant during the relaxation processes. The subscript e denotes 
the asymptotic limiting energies associated with equilibrium. 


Shock Wave Structure 

The previous examples, figures 1 through 6, rely on the "box" model, which 
is based on the assumption that the distributions have no spatial dependence. 

The Monte Carlo method, of course, has potential for much greater geneality. 

We can effectively introduce a spatial dependence into the distributions and 
study more complicated problems. To demonstrate this effect, results, of simula- 
tions for a normal shock wave structure are displayed in figure 8. For this 
example, the number of molecule in the sample size was not increased and, 
therefore, the curves are not exactly smooth. 

In this figure are displayed the translational and rotational temperatures 
(based on average energy) and density at seven distinct instants of time. As 
one might expect, the translational temperature develops an overshoot. As the 
rotational mode is excited the high translational temperature decreases and 
approaches an asymptotic steady value. 

This example is included to demonstrate that such simulations that involve 
both elastic and inelastic collisions are possible. More refined shock shapes 
than those displayed in figure 8 will require a considerable -increase in the 
number of molecules within the statistical sample and in computation time, 
thus, no attempt has been made to check the convergence of the solution. 


CONCLUSIONS 


The Monte Carlo simulation method described in this report, including the 
use of the Itikawa model for representing inelastic collision processes, is a 
viable scheme for studying translation-rotation interactions. The method can 
provide very useful qualitative and quantitative information on the relaxation 
processes associated with at least relatively simple topological systems (i.e., 
one-dimensional and quasi one-dimensional systems) . On the basis of experience 
gained here, it is not expected that the method will be considered viable at 
this time for more complex topological studies (i.e., three-dimensional flow 
simulations), because current and foreseeable computer resources appear insuf- 
ficient to allow economical processing of the increased sample size that will 
be required in such studies. 
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The method, however, is very useful in its present form to visualize 
fundamental gas kinetic behavior as demonstrated by the results presented in 
this paper. A review of the simulations given in this report shows t-he -fol- 
lowing results: 

1. Single step (Aj = ±2 for homonuclear molecules) transitions are the 
significant mechanisms of intermodal energy transfer rather than the multistep 
transitions (i.e,,. |Aj|>4 for homonuclear molecules). 

2. The coupling of translation-rotation transitions is the most effi- 
cient for low lying rotationally excited molecules and is least efficient for 
the highly rotationally excited molecules. 

3. The "relaxation time" required for molecules to reach an asymptotic 
steady-distribution in both the velocity and rotational states is dependent on 
the initial distributions. 

4. Relaxation occurs via a successive set of distributions that are not 
Maxwell-Boltzmann (nonlocal Maxwellian) . 

5. Initial rotational distributions with high rotational energy and 
that are far removed from satisfying equipartition lead to the appearance of 
"satellite peak" on the velocity distribution via a mechanism that is similar 
to the Stokes Raman effect accompanying the Rayleigh scattering. 

Subsequent studies should quantitatively compare characteristic relaxa- 
tion times found by the Monte Carlo methods with s imil ar times obtained 
experimentally. Of course, only qualitative comparisons are given here. 

The simulations reported in this paper certainly demonstrate that the 
method is viable for studying translation-rotation interaction processes and 
that some revisions are necessary in existing analytical methods. 
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APPENDIX A 


ESTIMATE OF SATELLITE PEAK POSITIONS IN THE VELOCITY DISTRIBUTION 


Certain types of rotational energy distributions can couple strongly 
through the collision processes to perturb the velocity distributions. In 
these cases, peaks -occur (on the velocity distributions) that are analogous 
to the Stokes and anti-Stokes rotational Raman effects which accompany Rayleigh 
scattering. The satellite peaks have been observed experimentally (refs. 16 
and 17). In one case (example 1) energy is "dumped" into a narrow band of the 
rotational energy levels; a pair of "satellite peaks" then appear around the 
maximum in the velocity distribution during subsequent, relaxation. In another 
case (example 2) a similar effect occurs when the velocity distribution has a 
peak at low velocities ("cold gas") and the rotational energy is peaked at the 
higher levels. Here, however, only one "satellite peak" appears in the 
resulting velocity distribution. The nature of these peaks is such that their 
position can be readily estimated without solving the Boltzmann equation. 

The relative velocity, g’, of a pair of molecules after a collision can 
be found from equation (8) and is given by 


g' = g /l - AE r /(l/2 Ug2) 


(Al) 


where AE r is the change in rotational energy during a collision, and p is 
the reduced mass. We assume that only one of the pair of colliding molecules 
transfers rotational energy during the interaction; the rotational transitions 
correspond to j -*■ j ± A, that is, 

AE r (j -»• j ± A) = E r (j ± A) - E r (j ) 

= A ■ (A ± (2j + 1) )k 0 r (A2) 

where k is Boltzmann's constant, and is the characteristic rotational 

temperature which corresponds to a single transition. We introduce the fol- 
lowing notation 


velocity corresponding to. the initial Maxwellian peak: c* 
velocity after rotational- translational interaction: c f 
reference velocity: Cq = /2kT e /m 

and, in addition, the approximations g ~ 2c* and g' ~ 2 c'; we then obtain 
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c» c* / A E r / f 1/2 u(2 c*) 2 ■ 

c 0 eg ]/ kT e j L 1/2 id cq 2 




| (A ± (2j +1)) 


Ti 


(A3) 


where we have used U = m/2 and (c */cg)^ = T^/Tg. We recall for homonuclear 
diatomic molecules, .such as molecular nitrogen, that multiples of | A j = 2 
rather than single transitions f A | =1 are allowed. In the analogy with 
Raman- scattering, the positive sign corresponds to "Stokes” and the negative 
sign to "anti-Stokes" effects. In the discussion above, example 1 displays 
both Stokes and anti-Stokes effects; example 2, however, shows only anti- 
Stokes effect. 
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APPENDIX B 


MODIFIED ROTATIONAL TRANSITION PROBABILITY 


In the Introduction, we listed several properties that the rotational 
transition probabilities must satisfy to ensure their proper behavior for the 
Monte Carlo simulation method: (1) probability must be conserved, (2) proba- 

bilities relating transitions to and from pairs of definite states must satisfy 
"detailed balancing," and (3) the probabilities, when used in the simulations, 
must yield the correct asymptotic behaviors of the distributions. The third 
property has been covered in the text (also, see refs. 4 and 5). The first 
property, that is. 


2 : 


t P(i, j -»- i', j gl = 1 


(Bl) 


is satisfied by Itikawa’s relations (refs. 5 or 11) as is also collision 
symmetry, given by 


(B2) 


To ensure the satisfaction of the second property listed above, we 
introduce the modified transition probability, P, given by 


P(i,j 


* * 

X J J > 


g) 


A (ca p fl 4 

(2i + 1) (2j + 1) ^ 


-*• i’ 


,j'; g) for (i, j) ^ ( 


i >3 ) 


(B3) 


Tliis relation also satisfies the "principle of detailed balancing" given by 

(2i + l)(2j + l)P(i,j i’ , j ' ; g) = (21* + 1) (2 j » + l)P(i»,j f -> i,j; g') (B4) 

.where A = A(g) and the symmetric function B(i,j; i'-, j’) are arbitrary rela- 
tions that have functional behaviors as indicated in the parentheses. Several 
example relations of B(i,j; i , ,j’) are 

B(i,j; i f ,j’) = [(2i< + l)(2j< + l)] a ‘ (B5) 

or 

B(i, j; i’,j’) = [ (2i + 1) (2 j + l)(2i’ + 1) (2 j f + l)] 1 / 2 (B6) 

The notation i< is used to designate the smaller value of either i or i’ . 
and similarly for j<. 
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Note that equation (B3) is obtained from equation (B4) , since g -»• g 1 
implies A(g) ~ A(g’). Itikawa's expression derived- in reference 5 is a 
sp.ecific -example of the more general probability relation that is displayed 
here (e.g. 5 Itikawa assumes A(g) = 1 and a = 1 in equation (B5)). It can 
be shown that the modified transition probability given by equation (B3) also 
satisfies "conservation of probability" and "detailed balancing." 

It is worthwhile to comment that equation (7) in the text, which is the 
selection rule for inelastic encounters, was based somewhat on heuristic 
arguments and yielded qualitatively satisfactory results. We expect, however, , 
that more accurate representations that will be based on more convincing physi- 
cal arguments, which will involve A(g) above, can be obtained for this 
equation. 


20 



APPENDIX C 


PROGRAM LISTINGS 


The entire program for the "Gas in an Imaginary Box" calculation is 
listed in this section. 

Program listings consist of sample control cards, correction cards, main 
program listings, and sample input-data cards. Many unused cards are still 
in the listings, but are marked by a comment symbol "c," "c*," etc. 
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reproducibility of the 

* ORIGINAL PAGE IS POOR 

* 

* SAMPLE CONTROL C.ARD5 

-J.- ' ’ 1 

* 

YOSHX, ”0500. ST0P9, X6363, YOSHIKAWA 

ACCOUNT, S”GKKY,T460<F>. 

AUDIT , ID= YOSHIKAWA. 

ATTA CH , TAP E 1 1 , FI L E702 , I D=Y0SHI K AW A, PW=STGKKY , MR= 1 , CY=2 . 

ATTACH, OLDPL, ITILIBS0URCF'»I0=Y0SHIKA-MA»MP*1. 

UPDATE f F, 

FTN, I »R=3»QPT=2»PL=l-00000. 

ATTACH, I MSL, TMSLL IB, ID= AMESLIB. 

LIBRARYUMSL) 

REQUEST, TAPE 9, *PE . 

REQUEST, TAP E10,*PF. 

L D SF T , MA P=X. 

LOAD, LGO » 

N0GO , MATN. 

RETURN, LGO* 

MAIN. 

CATALOG, TAPE9,DXXX903,TD=YGSHIKAWA,PVI=STGKKY,MR=1 ,RP=999 , CY=2. 
CATALOG , TAPE 10 ,0XX1 03, I D= YOSHIKAWA, PW=ST €KKY,MR=I,PP=999, C Y=2» 
AUDIT, ID=YnSHIKAWA. 

EXIT. 

C A T A LOG, “A PE9 ,0XXX903 , 1 C= YOSHIKAWA, PW=$TGKKY , MR= 1 , RP=999, CY=2 • 
CATAL0G,TAPE10,DXX103, TD= YOSHIKAWA, °VJ=STGKKY,MR=1,RP=999,CY=2. 
AUDIT, ID=YOSHIKAWA. 


* 

* PROGRAM LISTING 

if 

* 

DENT, CORRECT 
*7 INI *$.18 

C OTHER RANDOM GENERATION 

KRAN=50 

DO 50 JP-l , KRAN 
R=RANF (0) 

50 CONTINUE 

*PECK M0N T -5$ 
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REPRODUCIBILITY OP 'mu; 
ORIGINAL PAGE IS POOR 


PROGRAM MONTECC INPUT , QUTPdT ,T APE5=INPUT * TAP E6- OUTPUT »TAPE7»TAPE8» 
l T APE9-» TAPE 10 jTAPEII) 

COMMON /TIME/ T0 t TS,TF,TM t DTO*DTM,TN 
COMMON /CV/ MAX,MAX6,C1*RH01 
COMMON /CONST/ W, A, V0M*SF ,CO 
COMMON /RANDOM/ R 

COMMON /ANSWER/ MOM t t (4 ) , G ,DV0L , AI ♦ POA, ROT I 
COMMON /PART/ PC 5001) 

COMMON /CVl/ABCC(40,40,9) 

DATA ABCC/ 14400*0. / 


MONTE CARLO PROGRAM FOR GASES IN A BOX -REVISED BY K.K. Yoshfkawa 
JUMP=0 

1 CALL INPUT (JUMP) 

IFCJUMP.GT.O) GO TO 30 
10 CALL INITAL 
CALL MOMENT 
CALL OUTPUT CM ) 


30 

IF(TN.GE.TM) 

GO 

TO 

100 

50 

CALL JPAIRS 




100 

TM = TM+DTM 




200 

TFCTM.LT.TS) 

GO 

TO 

30 

209 

CALL MOMENT 




330 

IFCTM.LT. TO) 

GO 

TO 

30 


CALL OUTPUT ( M) 




TO=TO+DTO 




340 

IFCTM.LT. TF) 

GO 

TO 

30 


ENDFILE 10 
CALL EXIT 
STOP 
END 

*DECK IN P$$ 

SUBROUTINE INPUT (JUMP) 

COMMON /TIME/ TO , TS ,TF , TM- t DTO,DTM, TN 
COMMON /CV/ MAX,MAX6,Cl t RH01 
COMMON /CONST/ W, A,VOM*SF »C0 
COMMON /RANDOM/ R 

COMMON /ANSWER/ MOM ,T (4) , G, DVOL , AI , POA, R0T1 
COMMON /PART/ PC 5001) 

DIMENSION FCC40); FWC40) 

C 

C DATA INPUT FOR MONTE CARLO PROGRAM ALA BIRD 
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C 


DIMENSION HEDU8) 

REAO( 5,500 ) HED 

READX5 ,501 ) MAX-, RHOI, CO , W , A,VOM, MGM,DTM, TS, DTO, TF ,D,R 
10 MAX6=5*MAX 

AI=.25*W*D**2 
120 POA=SQRT (3.14159/A) 

200 CI=CO 

DVOL = 0 .5* (W/ (RH01*A) ) **2 
SF = RHOI* DVOL/ ( M AX*W ) 

WRITE? 6, 600) HED 

WRITE (6,601 ) M AX , RHO 1 ,.DTM»TS, DTO *TF , CO ,W* A , VOM, SF ,D 

500 FORMAT ( 18A4) 

501 FORMAT { 14/ 10X5E10.4, I3/f 4E10.4)) 

600 FORMAT (1H1 18A4) 

601 FORMAT?* MAX=*I 15//5OX*RH01=*G15.6 DTM=*G15 .6, *;TS=*G15.6, 

1 * DT0=*G15.6//* TF=*G15.6,* C0=*G15.6,* W=*G15. 6,*; A=*G15„6, 

2*;VOM=*G15,6,* SF=*G15.6//* D=*G15.6) 

IF (MAX.LT.O) GO TO 300 
RETURN 

30 0 READ? 8) MOM , T, G, DVOL ,W, A , VOM , S F,CO, TO, TS ,TF ,TM, DTO, DTM,TN, MAX, MAX6 
1 ,C 1, RHOI , P , A I , R0T1 , POA , A3CC ,R 
REWIND 8 

NTMO=? TM-DTM) /DT0+0.5 

100 READ(T) TAUBAR,NTM,MAX,C0C2,E0E,FC,FW 
WRITE ( 101 TAUBAR,NTM,MAX,C0C2, EOE, FC, FW 
IF(NTM.LT.NTMO) GO TO 100 
TF=TF+100.*DTM 
JUMP=1 
RETURN 
END 

*DECK T N I $ $ 

SUBROUTINE IMITAL 

COMMON /TIME/ TO , TS ,T F, TM , DTO, DTM,TN 
COMMON /CV/ MAX,MAX6,C1,RH01 
COMMON /CONST/ W, A, VOM , SF ,C0 
COMMON /RANDOM/ R 

COMMON /ANSWER/ MOM ,T(4) , G, DVOL ,AI , POA, R0T1 
COMMON /PART/ P( 50011 
COMMON/INI X2/ ERO 

INITIAL VALUE SUBROUTINE FOR EQUILIBRIUM BOX 
OTHER RANDOM NUMBER GENERATION 
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C ' KRAN=1 0 50 OR 100 ETC.' 

• DO 50 JR=]>KRAN 
R= RANF (0) 

50 CONTINUE 

ERO-O . 5*W*C 1 1 

J0=S0RT ( I .+*9940 1 E16*ER0 } / 2 • 

CM=0. 707106 B*C1 ( l/SQRT(2))*Cl 

ISEED=4111U111 
TN=0. 

P(1)=0. ‘ 

DO 100 I M= 1 » 4 
100 T( JH) = 0. 

CALCULATION OF INITIAL VALUES 
1=0 

130 pci)*p(n+i, 

DO 200 L=2.»4 
CALL RANDimX, TYtR) . 

I X=T Y 
IL=T+L 

200 CALL GAURND(IX,CH,0.,P(TU) 

200 PUL l=CM*GGNOFU'SEED) 

ROTATIONAL ENERGY AT UPSTREAM CONDITIONS 
CALL X2DTSTCIR0T) 

P(I+6)=IR0T 

400 P(I+5)=S0RT(P{ I+2)**2+P( I+3)**2+P{ I+4)**2) 

1 = 1+5 

TF{ I .GT.MAX6) GO TO 500 
GO TO 130 

READ INITIAL VALUES FROM THE FILE 11 
R'EADl 11) P 
500 CONTINUE 

DTO=nTn*D T H*. 99999 
TS— TS^DTM# . 99999 
TF=TF*D7M*.9999 
TO=TS+DTO 
TM = DTM 
RETURN 
END 

*DECK JPA$$ 

SUBROUTINE JPATRS 

COMMON /TIME/ TO , IS ,TF, TM ,DTO,DTM,TN 
COMMON /CONST/ w, A , VDM t SF ,C0 
COMMON /RANDOM/ R 

COMMON /ANSWFR/ MOM , T ( 4) ,G»DVOL , AT »P0A,R0T1 
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COMMON /PART/ P 15001) 

C 

C COLLISION PAIR SELECTION FOP. THE INVERSE P.Q.WE-R POTENTIAL '»'i~TT CIES 

C 

VP=.0 

10 VM-P { 5 ) 

NC=P(1 ) 

20 DO 40 I = 2 5 N C 
J=5*I 

TF(P{j},Gf.VM} GO TO 35 
IF(P(J l.GT.VR) VR^PCJ) 

GO TO 40 
35 VR— VM 
VM=P { J ) 

40 CONTINUE 
50 GM=VM+VR 
60 R=RANFCO) 

80 J1=P(1)*P.+ 1 
90 P, = RANF { 0 ) 

1X0' J2=P(1)*R+1 
120 IF ( J.l . EQ. J 2 ) GO TO 90 
J1=5*J 1-3 
J2=5*J 2-3 

150 G=SQRT((PC JU-PC J2) )**2+CP<Jl+l ) 

1-PCJ2+1) }**2+CP( Jl+2)-PC J2+2) ) **2) 

C PAIR SELECTION RULE: IF FCG/GM) > R, TAKE A PAIR 

C FOR GENERAL CASE USE CC1 THROUGH CC3 

C F { G/GM )- ( G/GM ) 

R=P.AN p (0 } 

FGBAR=CG/GM 1**0. 4310 
KELST=0 

C KELST= 1 FOR MONATOMIC GAS 

IFIR.LT.FGBAR) GO TO 170 
CC1 CELST^O .21 5 FOR EXAMPLE' 

CC2 FGEL= { G/GM ) **CEL$T 
CC3 I F{ R-. GT. FGEL ) GO- TO 60 

CC* THIS PROGRAM IS SET UP FOR MAXWELLIAN MODEL CASE (FGEL=l«0 ) 

KELST^l 

C END OF PAIR SELECTION RULE 

170 CALL CRASH(PCJi) , PC Jl+l) , PC Jl+2), Pt J 1+4) ,P(iJ2) , 
lPCJ2+l) t P(J2+2) »P (J2+41 f KGP.f KELST) 

TFCKGP.GE.l ) GO r n 60 

P(J1+3) = SQRT{P( Jl)**2+P( Jl+l)^2+P(JI+2 1 )^?) 

PIJ2+3 SORT C P { J2)**2+PU2+1)**2+P{ J2+2)**2) 

260 DTM=2 ,*OVOL/C PC 1)**2*A*G*SF ) 
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270 TN = TN + OTN 
280 TF(TM.LT.TM) GO TD 60 
290 RETURN 
END 

❖DECK CRA$i 

SUBROUTINE CRASH { U1 , VI ,WI , Il,U2fV2«W2f 12 KGP, KELST ) 
COMMON /CONST/ W* At VOM, SF,CO 
COMMON /RANDOM/ R 

COMMON /ANSWER / MOM , T{ 4) , G ,OVOL » AI , PDA ,R0T1 
CQMM0N/PRB0UT/PSUM(400) ,RRRR,.JIFIN, J2FIN • 

REAL II, T2 

C DIATOMIC PARTICLE COLLISION TRAJECTORIES 

C *EL5TIC COLLISION ; GO TO 700 

IE(KELST.EO.I) GO TO 700 
G2=0*G 

IF(T1.GT.38..0R.I2,gt.38.) GO TO 7 
ETP=0.25*W*G2 

C EXCLUSION OF NO ENERGY TRANSITION AT LOW KINETIC ENERGY 

GO to I 
7 KGP= 1 
RETURN 
1 P-RANFCO) 

EPS= 6 • 2831 8*R 
KG P=0 

ERl=4,0241E-16*Il*ni+l> 

ER 2=4. 0241E— 16* I2*( 12+1) 

RM1 =S CRT (2-*AT*ERl )' 

PM 2 =SQRT{2.*AI*ER2 ) 

10 R= RANF (0 ) 

50 V02= ( VOM**?)*R 
VO=SQRT( VO 2) 

60 CONTINUE 

IF (M0M.LT.100) GO TP 100 
R2=V02*A/3. 14159 
GO TO 110 

100 B2 = Vn?*A/G**(4./ M CM) 

ROTATOR MODEL 

TRANSITION PR0BA3ILTTY FOR ROTATOR l 
TRANSITION PROBABILITY FOR ROTATOR 2 
110 CONTINUE 

B8=S0RT I 82 1 
LI 0=1 1 + 0. 1 


ma-EODTJcroaire of ™ 

ORIGINAL PAGE ® 
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L 20= 12+0.1 
12 0 R=RANF (0 ) 

C SEE PROB. STATEMENT NO. 600 

•RRRR= R- 

CALL LINKMCIL10,L20, BS* ETR) 

J J 1= J 1FI N 
JJ2=J2PTN 

130 Jl=2* J Jl— 2+MOD (L 1 0*2 ) 
J2=2*JJ2-2+MOD(L20,2) 

# # # # # % 5}: ^ A s}c ■%. # # # iff. $ * * * ^ # s£ If % 

ER1P=4.0241E-16*J 1*{J1+1) 
ER2P=4.0241E— 16* J 2*1 J2+1 1 
200 GP2= (G2-4. *(£RlP + ER2P-{ ER1+ER2 ) )/W) 

IF (GP2 »LT .0.0} GO TO 150 
GO TO 160 

150 IF (KGP.GE.l) RETURN 
KGP=KGP+1 
GO TO 10 
160 GP=SQRT(GP2) 

RM1P=SQRT(2.*AI*FR1P) 
RM2P=SQRTI2'.*AI*ER2P) 

205 BP=( G*8B +2 . *( RM1+RM2-RM 1P-RM2P } /W) /GP 
IFI8P.LT. 0.0) GO TO 10 
GO TO 210 
TOO G2=G*G 
701 R=RANF (0) 

EPS=6.28?13*R 
710 R=RANF 10 ) 

750 V02=I V0M**2)*R 
VO=SQRT{ V02) 

760 CONTINUE 

IF ( MOM. LT .100) GO TO 800 
B2=V02*A/3. 14159 
GO TO 310 

800 B2=V02*A/G**{4. /MOM) 

810 CONTINUE 

8B=SQRT<B2) 

860 GP=G 
905 BP=BB 

J1=U+.001 
J2=I 2+ .001 

210 IF IMOM.LT.IOO) GO TO 220 
V0P=R p *P04 
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IFtVOP.GT.. *39995 } VOP=. 99995 
GO T 0 230 

220 VOP=RP*GP**I2,./MOM) /SQRT(A) 

230 XT=SCATER(Vn,HOH ) 

XI P= SC ATFR ( VQP ,MOM) 

250 XI8=.5*(XI+XIP ) 

260 GX=U2-U1 
GV=V2-V1 
GZ=W2-W1 
300 CE=COS (EPS ) 

SE=SINIEPS) 

CX=COS (X 18 J 
SX=SI N (X I B) 

GPG=GP/G 

350 RTG=SQRT(G2-GX**2) 

GPX=GPG*(GX*CX-RTG*SX*CE) 

G PY= G P G* ( G Y *r X+ ( G X *GY*C E+G#GZ*SE ) /R TG* SX ) 

GPZ=GPG*(GZ*CX + IGX*G7*CE-G*GY*SE >/RTG*SX) 

400 U2=» 5* {U1+U2+G PX ) 

V2=.5*(V1+V2+GPY } 

W2=» 5# IWI+W2+GP7) 

U1=U2-GPX 
VI =V2-GPY 
VJ1-M2-GPZ 
460 I1=J1 
I2=J2 
RE TUP N 
END 

*DEC K LINKS* 

SUBROUTINE- L I NKMC (1 1 » T 2 , BB ,ETR ) 

DIMENSION PWAVEI 20*20) »LABCSA( 16) 

COMMON /CM l / PWAVE,FKIN,L10,L20,BIMP*NMAX, NPRINT , L 1PAR ,U2P AR 
COMMON /CM 2 / JEL 1 *JEL2» LA8C SA* LLMAX*I SELCT 
COMMON /CM3/ VA» VR* IPFT 1* TPRT2* I PRT3* IPRT4 
COMMON /CMVR1/ VC *VALPHA ,VC6* BBI MP, EEEE 

C** INPUT ** SELEC T ION FOR OUTPUT (0 FOR PWAVE/1 FOR P) 

I S ELCT -I 

C** INPUT ** MIN PROBABILITY 

C** INPUT ** INITIAL ROTATIONAL STATES 
LIO=! 1 
L20=I2 

C** INPUT ** RELATIVE KINETIC ENERGY UN EV) EKIN= ETR/1.602E-12 
EKTN=ETP*. 624219 7E+12 

C** INPUT ** IMPACT PARAMETER CIN ANGSTROM) 


REPRODUCIBILITY OF THE 
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BI'MP=RB*1. FOB 

C** INPUT ** MAX NO. OF TERMS IN EXP 
NMAX= 30 

C** TN PUT INDEX FOR PRINT OUT 
I PR T 1= l 
IPRT2= 1 
JPRT3= NMAX+1 
I'PRT4= 0 

C** INPUT '** POTENTIAL PARAMETERS FOR SPHERICAL PART 
C** VC R ) = VC*EXP(-VAL D HA*R)-VC6/R**6 

C** V IN EV » R TN ANGSTROM ' 

VO3440. 

VALPHA=3 .160 
VC 6= 73 .40 

C** INPUT ** POTENTIAL PARAMETERS FCR NON -SPHERICAL PART 
VA- . 2 
VB-.2 
JMAX= 20 
5010 CONTINUE 

DO 5016 I=1,-JMAX 
00 5016 J=1,JMAX 
PWAVEC I» J)= 0.0 
5016 CONTINUE 

c******** nprint 

NPRI NT= NM AX 

IFCEKIN. LE.O.OOl) EKIN= 0.001 
IF CEKI N.GE. 0. 5> EKIN= 0.5 
CALL PROP 
5999 CONTINUE 
RETURN 
END 

*DECK PRO$$ 

SUBROUTINE PROS 

COMMON/PRBOUT/ PSUM (400) ,RRRR,-J.1'FTN,J2FIN 
COMMON /M V 1 / AMATRX(?0,20,9) 

COMMON /MV2/ VBBf VAA, RROT » ETOT ,BRC , WALP 
COMMON /CM VR1/ VC ,V ALPHA ,VC€, BT MP f .EEEE 
COMMON /MV 3 / NC01JNT 

COMMON /CM1 / PWAVEfEKINf LL0 t L20» BBBB,NMAX-,,NPRTNT*L 1PAR ,L2PAR 
COMMON /CM 2/ JELI , JEL2 ,LABCSA, LLMAX,. I SELCT 
COMMON /CM 3/ VA, VB, TPRT1, TPRT2. TPRT3 , IPRT4- 
DIMENSION PWAVEC (20,20)" 

DIMENSION POOD (20,20), PE VN ( 20 ,20) , PWAV.F ( 20*. 20 ) 

DIMENSION AKSUMO ( 20, 20 ) , AKSUM V( 20, .20) 


30 



DIMENSION L ABCS A { 16 ) 

TCLOCK= 0 
BI MP=8BRB 
JMAX=?0 

LMAX= 2*JMAX-2 
C********* R C DUC ED M^SS 
RMA$$= 14.02 

****** ROTATIONAL CONSTANT 
3R0' r * 0.2512 F~3 

VVAL P= 0.045723*VAL D HA/ SQRT< RMASS) 

VAA— 0 .447 2 1 36*V A 
VRB= 0.2*VR*0. 6298283 

ETOT= EKIN + FLOAT CL 10*(L 10+1} +L20* (L20+1 ) ) *8R0T 
E8= ET C T /B ROT 
FBI* SORTTEB) 

LLMAX= I NT f FBI )+l 
£***** PRINT qf)l-l 
11 CONTINUE 

TFfLLWAX.GT.lMAXJ LLMAX=LNAX 
DO 13 T= 1 i JMAX 
DO 18 J= l» JMAX 
PEVN(IfJ)* 0.0 
PODDCI ,J1= 0.0 
PWAVE(ItJ)* 0.0 
PWAV EC ( I » J) = 0.0 
AK SUMO ( I, J >=0.0 
AK.SUM1 ( T , J 3 = 0.0 
DO 16 M*l,9 
AMATRXU, J,M)« 0.0 
16 CONTTNUF 
18 CONTI NUF 
LL10= L 10+ 1 
LL20= L20+1 
JFL1 = (LL10+D/2 
JFL2= C LL 20+ 1 ) / 2 
LI PAP = MOD (L10»2)+l 
L2PAR= MOD ( L20 » 2) +1 
CL0= FLOAT! < 2*L 1 0 +1 > * ( 2*L 20+1 H 
PE 0= 0.0 

£********* J*J= i 
N= 1 

LLII* LL 10 
L L 2 1 = LL 20 
NCOUNT* 0 
PQ 50 K= 1 1 9 
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CALL VMATRX (LL1T, LL2T t K) 

5° CON T imjF 

LL1FMM= MAXOCl L10-?,L1PAR) 

LLl-FM-X- LL 1:0 +'2 

LL2FMN= f'AXOI LL20-2 ,L2PAR) 

LL2FMX= l L 20+ 2 
K— 1 

J J T 1 = ( L L 1 1 + 1 ) /•?. 

JJI? = (LL2T + 11/? 

DC! 64 1=1,3 
LL2F= LI 20-4+2*1 
DO 63 J= 1 » 3 
LL 1F= LL 1 0— 4+2* J 

TFRL2F.LT. 1.0P.LL1F.LT.1) GO TO 62 
JJF1= (LL1F+D/2 
J JF2 - (LL2F + U/2 

AKSUMO(JJFl f JJF2)= AMATRX( JJI1» JJI.?,K) 

62 K=K+1 

63 CONTINUE 

64 COMTTNUF 
TOTPW= 0.0 

DO 74 LLIF=LL1FMM,LL1 C MX,2 
00 73 LL2F=LL2FMN t LL 2FMX t 2 
J JF 1 = RLIF + 1 ) /? 

J J F?= (LL2F + D/2 

POOD { J JFl » J JF2 ) = AKSUMO ( JJF1, JJF2) 

P F VN ( J JF 1 1 J JF2 I — 0.0 

IF (l L1F.F0.LL1I. AND.LL2F.FQ.LL2T ) PFVNUJF1 ? JJF2)=1.0 
PWAVEI JJF1, JJF?)= PODO(JJFl,JJ F2 )**2*4.0+PEVN( J J>= l, J JF2) ** 2 
o=PWAVEf JJF1, JJF2 ) 

LO 1 = L L l p -1 
L02=LL?=-1 

(■/***** PRINT 902-1 
7? CONTINU- 

TQtdvj=totPW+P 

73 CONTINUE 

74 CONTINUE 
G*#*** ppjmt 904-1 

75 CONTI N Ur 

IF(NMAX.EO.l) GO TH ^?00 
N= 2 

C*^*«^«***100 
100 C OM~T NUE 
NCOIJNT* 0 
M2= 
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L L 21 MM= MA XO ( LL 20 - N2 , L2 PAR )• 

LL2 I MX= LL20+M? 

! F t LL OF .L AX) LL?TMX= LL MAX 
DO 140 U 2T»LL2TMH,ll.2I«X*? 

JJI2= (LL2I+1J/2 

cp, p= v^~(eb- c L oat(lL’T-1 )**?. \ 

llOAX’s INT(=3?)M 

LLL T = LL10+N2 

JJIl={ LlU + U/2 

IF (LLMAX1.GE.L.MAX) LLMA-Xl= LMA X 
TFCLL1 1.OT.LLMAXl } GO TO 130 
rn 114 K=5, a, 3 
CALL V MATRX (LL1T » LL 2 1 » K ) 

114 CONTJ NOE 

00 IL G K=2, 5,3 

TALL V^ATRXILLIT, LL?T,K) 

■119 CONTIMIJ- 

1= (LL 1I-2.LT. 1 ) OQ in 125 
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AMATFX (JJTl , JJ T2 ,4) = AMATRXUJIi-l, JJT2,A) 

K=4 

TFILL2I .PO .I.L20-N ?.np . LL2 1 .EQ .LL 20+N2 ) CALL VMA TRX ( LL1I ,LL2I,K5 
AMA.TPXUJT1 , JJI2 , 7}= AMA-RXIJJ 1 1-1 , J J 1 2 + 1 , 3 ) 

K- 7 

IF ( LL? T * c 0 .LL 20 +N 2 - 2 .OR .LL 2 T.E 0 .LL 20 +N 2 ) CALL VMATRXI LL 1 I ,LL 2 I ,K) 
TFILL 2 T- 2 .LT. 1 ) f -0 TP 130 

a«AT«>x( JJT 1, JJI?. l)= AMA T RXf JJT 1-1 , J J 1 2-1 »9) 

K= l 


T F ( LL2 I .BO. LL20-N2+2 . OR . LL2 I . HQ. LL2G-N2 ) CALL VMATR X ( LL 1 I » LL2I ,K 1 
125 CONTINUE 


TF (LL2 I-" 1 . L T .l ) GO TO 130 

AMA T PXf JJI1 , JJI2.2) = AMATRXl JJT1,JJI2-1 ,8) 

K= 2 

TF(LLOT.FQ.LL?0-N?) CALL VMATRXILL1 1,LL2I ,k> 
1.30 CON’INUF 


LLLT=LL10-N? 

JJI1= (LLlT+l)/? 

IF (LL II. LT.L IP ftp I GO tq 149 

rn 134 K=5 ,8,3 

CALL VNATRXULH * LL2I ,K) 

1?4 CONTINUE 

DO 139 K= L , 7 » ? 

CALL V M A7R X { LL1 1 , LL2.I »K) 

1 29 CON” T NUF 

AMft"PX( JJT1,JJ I? , 6)= AMATpxC 12,4) 
K= A 
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IF (LL2 T.EQ.LL20-N2. OR. LL 21.50. LL20+N2) CALL VMATRX(LLIT,LL2I » K) 
AMA T nX(JJUtJJI?»Si» AMATRX (JJI1+1 ,JJI2'+l,l ) 

_ I E LLL2 L. EQ...LL20 +N-2-? . HP .L L2T . EG . LL 20+N 2 ) CAL L VMATp X{ LL 1I,LL 21 ,K) 
IF (LL 21-2. L T. 1 } GO TO 149 

AM ATRX (J JT1, J J 12, 21= AM ATRX < JJI 1,JJ1 2-1, 8) 

K= 2 

IF(LL2I.FQ. LL 20-N 2} CALL VMA T RX{ LL1I ,11.2 I ,K1 
A?^ATRX f J JT 1 ,JJT2,3)~ AMATRXCJJI 1 + 1, JJ T2-l,7) 

K — 

T F ( L L 2T . FO.i L20-N2+2.0R. I.L2I . FQ.LL20-N21 CALL VMA TP X ( LL1 1 , LL2 I , K ) 
149 CONTINUE 
A**#*# i[50 

LL1 I MN= MAXO (LL 10— N 2+2, L IP AR ) 

LL17MX= LL10+N2-? 

TF(LL1TMX.GE.LLMAX1 LL1IMX= LLH4X 
DO 199 l. LI 1 =L L 1 1 M Ht LL IT MX , 2 
JJ T 1= (LL1T + D/2 
FI32= 3QP“{FB- FLQA’i'ILLlT-l} *$2 ) 

LLMAX1=TMT ( 9R2 ) + l 
LL2I= LL20+N2 
JJ 12- (LL2T+1)/? 

IF (LL^AXl .GF.LMAX ) LLMAX 1= LHAX 
JFLLL2I ,C T . LLMAX1 1 GO to 180 
DO 164 K=5,9 

I F ( LL1 1 • FQ» LL 10-+N2— 2 . AND. K . EQ .6 ) GO TO 163 
CALL V V ATRX (LL1I ,LL2I»K) 

GO TO 164 

163 rOL’-TNUE ’ : ’ 

AMA"PX(JJT 1,JJT2,61= AMAT.RXI JJI 1+1 , J JI2 ,4) 

164 CONI IN Ur 

IF (LL2I-2.LT. 1) GO TO ITS 

AMATRXf JJI l,JJI2,2)= AP.ATRX (JJI1,JJT2-1,8) 

AM'ATRX ( JJT1 ,JJT? ,3)= AM ATRX C J J 1 1 +1 , J J 1 2-1 ,7) 

IF(LL 1I-2.LT.1 ) GO TO 130 

AM ATRX (JJI1 ,JJI2.,1) = AMATRXC J J 1 1-1, J J J 2- 1,9) 

175 CONTMUF 

IF (LL 1. 1-2. L T.l ) GO TO -180 

AM ATRX (JJ!1,JJI2»4)= AM ATRX (JJ 1 1-1, JJ 12, 6) 

130 CONTINUE 

LL2 T = l L 20-N 2 
JJT2= (LL2I+U/2 
IF(LL2I.L T .L2PAR1 GO TO 199 
pn 1R4 k=5 , 6 

IF (LL1T.E0.LL10+N2-2.AND.K .E0.6 ) GO TO 183 
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CALL VMATRX (LL1T»LL?I*K) 

GO TO 184 

133 CONPINUF 

AMATRXtJJUtJJI?, 61= 4MATRXI J J 1 1+1, J J I?, 41 

134 CONTINUE 

00 189 K=l»3 

CALL VMATRX ( LL 1 T » LL2I ♦ K ) 

18 c CONTINUE 

AMATRX (JJIl»JJI2» 31= AMATPXC JJI1 , J JI 2+1 ,2) 
AMATRX { J JT 1 » JJ I2»9) = AMATRX (J J 1 1 + 1* JJ I 2+1, 11 
IF (LL1 T-2. LT» l ) GC T 0 1^9 

AMATp X ( J JT 1 1 J J 12 *4)= AMATBXt JJU-lt JJT2.6) 
AMATRX (JJIliJJT? ,71= AMATRX (JJ 1 1-1, J J 1 2+1, 3 ) 

C *i< 

199 CONTINUE 

LL1FMX=LL10+2*N 


LL 2FMX= LL.?0+2*N 

LLl FMN= LL 10— 2*N 

I C (LL1FV1N.LT.L1PAR) GO TO 201 
LL ll=LLl c MN+2 

200 CONTINUE 
LL2FMN=LL20-2*N 
IF(LL2FMN,LT.L2PAR) GO TO 202 
LL22= LL2FMN+2 

GO TO 205 

201 CONTINUE 
LL1FMN= LI PAR 
LL 11= L1PAP 
GO TO 200 

202 CONTINUE 
LL2FMN=L2PAR 
LL2?= L2PAR 

2 05 CONTINUE 

LL22F= LL2FMX-2 
IFCLL2FMX.LE. UMAX) GO TO 206 
LL2FMX=LLMAX 
LL2?F= LLMAX 
206 CONTINUE 

DO 249 L2=LL22,LL22F,2 

LL 11F= LL1EMX-2 

E8?= 50RTIE8- FLO A” (L2— 1 )#*?) 

LLMAX 1= INT(FR2)+1 
TF (LLMAXl.GE.LMAX 1 LLMAX 1= LMAX 
IF(LLlFMX.GE.LtMAXl) LL11F= LLMAX1 
DO 248 Ll= LL11,LLUF,2 . 
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Kl= 1 
K2= 0 
L2 ?=[.?-?. 

Ill- L 1-2 

IFCL2.I.F.L2PAR) L 27- L2PA* 

T p (Ll .LF.LIPA 1 ? ) L 1 1- L1PAR 
IF(L?.LE,L2PA2) K 1=4 
IFU 1.LE.LIPAR) K2= 1 
K= K1 

Jl= { L 1+1 ) / 2 
J2= (L?+l)/2 
AK= AKSUMOt J1»J2) 

L 2 D 2=L 2 + 2 

nn 219 L L2= L22,L2P2 t 2 
K = K+K? 

L 102=1 1+2 

DO PIS LL1=L11,L1P2,2 
JJ1= ILL1+U/2 
JJ?= (Li 2+1) /? 

AKSl'«l ( JJ1 1 JJ2 )= AK*4MATRX{J1,J2»K)+AKSUM1( JJl,JJ2) 

K= K + l 

21 8 CONTINUE 

219 CON TT NUt 

C********* OPT ’IT OF A'-jA^RX 
229 CONTINUE 
PAP- FOVTJ N0 C 
249 cnr-riMijo 
C ##* 2 5 0 

LLl p M X= LL10+2*N 

t otpc= o.o 

TOT PW = 0.0 

no 299 Ll? p =LL2FMNtLL2FWX,2 
LL 11 F= LllFMX 

F92= S0P.T(‘F3- FL n AT (LL2F— 1)**2) 

LLMA X 1= TN r (F3?)+l 
IF(LLMAXl.Oe.LMAX) LLMAXl= LMA X 
TF( LL1FMX.GE.LLM AX1 ) LL11F= LLMAXl 
nn 293 Ul p =H.tFMN,LlllF,2 
J J F 1 = (LL1F + D/2 
JJF?= (LL2F + D/2 

AKSHMK JJ C 1 , JJF2)= AK5UMK JJF1 ,JJF2J / FLOAT f N) 

T c (MOO CM, 2 ) .NE.O ) Gn TO 29 4 
AKSOMUJJPl ,JJF2)= -4.0*AKSUM1 (JJFl*JJF2 ) 

D FVM ( JJ P lt J J r 2 ) — PEVNCJJF1 T JJF2)+ AKSUMK J JF1 * JJF2) 

PUJAV F { JJF1 , JJF_?) = POOD (JJFI, JJF? )**2*4.Q+PFVN1 JJFl, JJF21**2 
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roNTiMur 

D POO( JJ C 1 , JJF2)= 3n r>D{ JJF1,JJF2)+ AKSU m ] (JJF1,JJ = 2» 

PWAVt ( J.l-l t JJ r 2)= PDnoiJJFl, JJF2)**2*4.0+P = VN( JJF 1, JJF2)**2 

295 CONTI NUF 

T= "VAVF* JJF1 , JJF2) 

LO 1= LL1F-1 
LO 2= LL2F-1 

O********^ DpJNJJ f)f p WAV? 

296 CON T INUE 
C*#*** PFIM“ 902-2 

297 CON T IN U- r 
TOTPW= TOTPW+P 

AK d )N0 ( J JF I , J JF 2 ) = AK5 U J JF I , J J F2. ) 

AKSU M 1(JJF1,JJF2I= 0.0 
2^8 CONTINUE 

299 CON-INUF 
PE1=PWAVE( JEL1, JEL2) 

n E 10= ARE( (9F1-PF0I/PE1 ) 

C *■**#* PR I N‘ r 904-2 
7 904 CONTI N!JF 

TP (AB^(Tn-PW-l.O) .f.O.lF-3. AMD.PF10.LT. 0.12-3) GO TO 300 
I F(N. EQ.NMAX) GO _r G 300 
N= M+l 

ocq = OF] 

GO 100 

300 CONTINUE 

£$«$$$$$$$ P[_ A c 71 C 

C********* FINAL PRINT 

I p ( I PRT4. FQ. 1) GO TO 325 

MDRI.r!T=0 
p.M AX= 0.1F-5 

301 C OM”T MLF 
C***** PP tnt 901-2 

ppjf.jT 910-1 
!«N= LL2FMN 
IMX=LL2PMX 

T F [MOO ( LL2-FMX— L? D AR, 2 ) . N F.O } IHX= LL2FMX-1 
CP 719^T = IMN,imx, 2 
LL?== r-’X + IMN-T 
J2= (LL2F+1)/? 

I L 11 F= L H F N X 

EB2= SPPTIEB- FLOAT! LL2F-1 )**2) 

LLMAX1= ? N T ( FB2) + l 
IFC1LMAXI.G5.LMAX) LLMA Xl= LMAX 
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JF(LLlFMX.Gc.LLMAXl ) LL11F=LLMAX1 
J1MX= (L L 1 IF+T ) /2- 
IF<J1MX.GT. 16) J1MX=16 
-L02-tt2F-l 
C**#** PRINT gti-i 

T P { m p p IMT. EQ * 1 ) GO TO 319 
LMIN2* MTNOC L20,LQ21 
C L 10= c LOAT ( 2*LM,IN2+I) 

J1MN={ LI PAR-H) /? 

nn 3i8 ji=j inn, jinx 

LO 1= 2 *U. 1-3+L IP AR 
LMIN1= MTNOM10,L01) 

CL 1= CL10* -L n AT{ 2*L f ^I N 1+1 ) 

P= PWAVEfJl , J2) 

PCOPCT= CL1/CL0*P 
TOTPC=TOTPC +PCORCT 
PWAVFCCJL, J2)= PCCRCT 

IF( J1.E0.JFU.AND.J2.EQ.JEL2) GO TO 318 
IF(ISFLCT.EO.O) pdut=p 
IF CIS cL<“T.E 0.1) PGUT= pcorct 

318 CONTINUE , 

319 CONTINUE . 

L ABC 9 A ( 1 ) = LIPAR-i 
DO 3 20 T = 2, 16 

LABCSA (T)=L,APCS A(1)+(I-1 )*2 
32 0 CONTINUE 
C***** PRINT 912-1 
C IF (sprint. EQ. 1 ) GO to 9906 

C***3s***** PR I NT OF PCORCT WITH PFLASTIC MODI FIFO 
325 CONTINUE 

PC= 1.0-TOTPC+PWAVECC JEL1 , JEL2 > ‘ 

PWAVECIJFLIt JEL2)-= PC ' - 

IF(ISFLCT.EC.O) PSUM(1)=PWAVE{ JEL1, JEL2) 

I c ( ISFLC T .EO.l) PSUMU»=PWAVEC(JEL1,JEL? ) 

STATISTIC MODI c I CAT ION OF ELASTIC COLLISION 
RRRP : P ( M) = P ( N)— P$U M ( 1 ) / ( 1 .—PS UN ( 1 ) ) PBAR VS RP.RR 

600 RPRP= { l.O-PSUM ( 1) )*RRRR 

IF ( PRRR . GT .PSUM ( l ) ) GO TO 1905 
J1FIN=JFL1 
J2 P IN=JEL2 
GO TO 2001 
1905 CONTINUE 
NN* 1 
N=1 

1«10 GON T I NUE 
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J1MN= JFL1-N 

TFt JIMN.LT.l) JX MN=1 

J1MX= JELl+N 

IF( J1MX.GT, JMAX) J1MX= JMAX 
K=0 

J2=JEL2-M 

IF { J2 .LT.l ) GO TO 1942 
192 0 CONTINUE 

00 1940 Jl= JlHNtJlMX 
NN=MN+1 

IF( TSELC7.EQ.1) PWAVE ( J1 , J2 )= PWAV EC ( J1 , J2 ) 
CALL SUMPIJ 1,J2,NN) 

C TFCRRRR.GT. PSUMINN) ) GO T 0 1940 

IFIPRRR.GT. (PSUM( NN)-PSUM(l) ) > GOTO 1940 
JlFIN^Jl 
J2FIN=J? - 
GO T 0 2001 
1940 CONTINUE 

IFIK.EQ. 1) GO TO 1950 
1942 K= 1 

J2= JEL2+N 

IF(J2.GT.JMAX) GO TO 1950 
GO TO 19 20 
1950 CONTI NUF 

J 2MN= JEL2-N+1 

IF (J2MN.LT . 1 ) J2MN= 1 

J2MX= JFL2+N-1 

IF { J2MX.GT . JMAX) J2MX= JMAX 

K= 0 

Jl= JFL1-M 

IF(Jl.LT.l) GO TO 199? 

1970 CONTINUE 

DO 1990 J2= J2MN,J2.MX 
NN= NN+1 

IF CTSELCT. EQ.l ) PWAVEI J1 , J2) = PWAVEC 1 J1 1 J2 )• . 
CALL SUMP( J1,J2»NN> 

C IF ( RRRR . LF • PSUM (NN ) ) GO TO 1935 

TFIRRRR.LE. (PSUM( NN)-PSUM(l) ) ) GOTO 1985 
IF(AB$U.O-PSUM(NN> J.LT.U0E-4) GO' TO 1985 
IFINN.GF.400) GO TO 1985 
GO TO 1990 
1985 CONTINUE 
J1FI N= J1 
J2FTN=J2 
GO TO 2001 


39 



1990 CONTINUE 

IF(K.EQ.l) CD TO 2000 
1992 K=1 

J-l-= JEL1+N 

I c ( J1 .GT.JM4X) GO TO 2000 
GO TO 1970 

2000 CON T I NUE 

N=N+1 

GO T P 1°10 

2001 CONTINUE 
1000 CONTINUE 

RFTlIRN 

END 

*DECK SUH$$ 

SUBROUTINE SUMPl Jl, J2»NNN) 

COMMON/PRROUT/ P SUM { 400 J » RRR R * J IF I N , J 2F TN 

COMMON /C. Ml/ PWAVE ,EKT N » 110 ,L20 ,.BB9B, NMAX , NPRINT , L IP AR , L2PAR 
DIMENSION P WAVE! 20* 20) 

N-NNM 

P?|JMCN|= PSUMIN-1 )+PWAVEf J1»J21 

RETURN 

END 

^DECK VM,AT*$ 

SUBROUTINE VMATRX (LL1 , LL2, K 1 
C#*#** PE VI SEP FQP ROOT 

COMMON /MV 1/ 4MATRXI 20 t 20,9) 

COMMON /MV2/ VBB ,VAA, BROT » ETOT , BRC, VVALP 
COMMON /CM VR 1 / VC ,VALPHA , VC6 , BI HP, E I J 
COMMON /MV3 7 N COUNT 
C ON MON /C VI / ABC C ( 40 f 40 , 9 ) 

L 1 J~ I.L1-1 
L2 J= LL2-1 

TPCK.GT.3) GO TO 1011 
L2I= L2J-2 

IFIL2I.LT.0) 1 GO TO 1510 
GO tq 1100 

1011 IFfK.GT.6) GO TO 1012 
L2T = L2J 

GO TO 1100 

1012 L2 T = 12J+2 
1100 CONTINUE 

IF{M0D(K,3) .NE.l ) GO TO 1111 
LI 1= L1J-2 
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IFU1T.LT. 0) GO TO 1510 
GO TO 1200 

1111 I F ( ‘TO 0(Kt3).NE«2) GO TO 1112 
LX 7= L1J 

GO TO 1200 

1112 111= L1J+2 
1200 CONTINUE 

ABC= ABCC(LLl,LL2,K) 

T-(ABC.NF.O.O) GO tq 1299 
CALCULATION OF VEFFI III ,L2 I / L1J Y L2J 1 
CC1= CG20C L IT tLIJ) 

CC 2 - CG20U.2I,L2J ) 

IFIK.LE.5) CSIGN= (-1.0)**{ L1J+L2J ) 

IFIK.GE.61 CSIGM= (-1.0)** IL1I+L2I 1 

C= PtnATf { 2*L1T+1)*( 2*L2I + 11*( 2*L 1J + 1 1 * ( 2*L2 J+ 1 1 1**0. 25*CSIGN 

B= VBB*CC1*CC2 

A=Q.O 

I = IL2I .NE.L2J ) GO TO 1215 
A= VAA*CC1/ SQR T ( FLOAT (2*L2T+1 1 ) 

IF (MOD TL 27, 2) .NF.01 A*-A 
1215 IFUir.NE.LU) GO TO 1219 

AA= VAA*C.C2/ SQ2T ( FLOAT 12* L1I +1 ) 1 
IF (NOD (L II, 21 .NF.O) AA=-AA 
A= A+AA 
1219 CONTINUE 

ABC= C B+A 1 
ABCCUL1,LL2,K)= ABC 
1299 C. ON T I N UF 

WI - BROT * FLOAT {L1I*{L1T + 1)+L2I*IL2I+11 ) 

WJ=RRCT* FLOAT (L 1J* ( L 1 J + l ) +L2 J *( L2J +1 1 ) 

WTJ= ABG(WI-WJ) 

E 1= FTOT-WI 
FJ = FTOT-WJ 

I F ( E I .LE .0.. O.OR . E J-LE . 0. 0 1 GO TO 1510 
EIJ= 0.5*{EI+FJ) 

CALL POOT(RC) 

RRC= 1 *0-{ 3 IMP/RC )**2+VC6/EIJ/RC**6 
El Jl = E7J*BRC 
I F { K .EQ. 5) GO TO 1500 

IFIL1I.EQ.L2J.AND.K.EQ.31 GO TO 1500 
IF (LIT .EQ.L2J.AND.K.EQ.71 GO TO 1500 
AIJ= ‘VVALP* SORTIE] JU/WIJ 

ew= eui/wij 

08 ALPH= EW/ATJ 
APAI= 1. 570796327 /A I J 
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F= EXP(-APAI) 

C A I J= 2.0*APAT*F/U.0-F*F) 

AAA= ABC^DBALP H*FAIJ 
149-8- CQNTI-NUE 
C*** 999 C-HFCK PRINT 

999 CONTINUE 

NC OUNT= NCOUNT+1 

1499 CONTINUE 

JJ1= (LL1+D/2 
JJ2= ( LL2+1 ) /? 

AM AIR X ( J J 1 1 J J2»K ) = AAA 
RETURN 

1500 CONTINUE 

AAA= ABC/VVALP* SQRT(EIJl) 

00 TO 1498 
1510 C ON T T \Mj c 
AAA= 0.0 
00 TO 1499 
FNP 

*D EC K -CG$S 

FUNCTION CG20( J1 , J?) 

C DOUBLE PRECISION FUNCTION CG2.0(J1,J21 

If (J2.E0.J1+2.OF .J2.EQ.J1-2) GO TO 8001 
IF { J2 . EQ. J 1 I 00 TO 8002 
C = 0.0 
on TP 8100 

8001 C0N T JNUF 

I F C J2 .E0.J1+2) J=J1 . 

I r f J2 .FQ.J 1-2) J=J? 

Xl= FL0A“(J+2)/ FLOAT C2«J+5) 

X 2= C L0 AT { J + l ) / =L04 T (2*J+3) 

X3= 1.0/ FLOAT (2* J+l ) 

C = RQR~ 1 1 « 5*X1*X2*X3 ) 

GO Tp 3099 

8002 CONTINUE 
J= J1 

Xl= FLOAT { J+l } / FLOAT ( ?.* J+3) 

X?= FLOAT ( J ) / FLOAT (?*J+1) 

X0= 1.0/ FLOA T { 2*0—1 1 
C= - S0RT(X1-X2^X3) 

80 n ? TF{M00(J.2) .ME.O) C= -C 
Binn C 020=0 
RFJURN 
FND 
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*PECK R00$$ 

SUBROUTINE ROOT! PC) 

C***** REVISED 3/26/74 

COMMON /CMVR17 VC ,VALPHA f VC6, 8T«P t EK IN 
RHIN= 1.12 

RRO= ALOGCVC/EKTHWALPHA 
I F (RP O.GE .4. 1 ) RRO = 4.1 
N=1 

3099 CONTINUE 
RR =RR 0. 

3100 CONTINUE 

RL= FKTN*BTMP**2/RR**2 
V= VC* FXP(-VALPHA*RR) 

Vl= -VAL PHA*V 
-If (VC6.FQ.Q.0) 00 TO 3101 
VP = VC S/PR**6 
V= V— vs 

V 1= Vl+6 . 0*VR/P.R 

3101 CONTINUE 

F= fV+RL-EK IN ) / ( 2 .0*RL /RR-V 11 
IF (ABSIF/RR) .LT.O .1E-5J GO TO 3199 
IF(N.GE.IOO) GO TO 3299 
RR= RR+F 

i‘F (BP .L T .RMIN) GO *0 3900 
N=N+1 

GO TO 3100 
3199 CON T I M UF 
RC= RP 
RETURN 

32 99 CONTINUE 

WRITE ( 6. 998 ) RR,F,RRO 

998 FORMAT { 1H0 //5X , 14HERR0R N GT 1 00, 3X, 3HRR=E13.5* 3X, 2E 13.5// ) 
RR=BPO 
GO TO 3199 
3900 CONTINUE 

P.RO= 0. 9*( PRO+RMI N.) 

GO TO 3099 
END 

*DECK OUT$$ 

SUBROUTINE OUTPUT (M) 

COMMON /TIME/ T0,TS,TF,TM,DTQ, D 7M f TN 
COMMON /CV/ MAX,MAX6,C1,RH01 
COMMON /CONST/ W , A , VOM , S F , CO 
COMMON /RANDOM/ P. 
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common /ANSWER/ HHM »T(4),G»QV0L.»AT, POA, ROT! 

COMMON '/PART/ o ( 5 00 1 > 

COMMON /C V l / AB CC. ( 40. t 40 ,9) 

-P-T-MF-NS'I ON FC {40 ) , FW (40 ) 

C 

MT-MM TM-DTM )/DTO + .5 
T&llBAr.=0.1^NTM^DTf3/ DIM 
f.nC»TI2»/{ T a)*Cl J*. 886227 
COC2=T n )/ (T( 1)*1.5*C1*C1 } 

EQE= . 8 0402 E-l 5 *“ ( 4 ) / (T { 1 ) *w*C 1 *C l ) 

WRITE (6,600) CO C. ,C0C2,E0E , T AUBAR 
DO 100 1=1,40 
C C(T ) =0 
100 B'J ( I ) = 0 

€0E=1 .0/Cl 
PO 200 !=5,MAX6,5 

jc=p{ n^coc^io. +i 
jw=p(i+u+i.oi 

I F ( J'C .GT.40 ) JC=40 
I P { J W . GT . 4 0 ) JW=40 
FC( JC)=FC( JCI+l. 

FW( JW)=FWC JW) +1. 

200 CnNTTMlJE 

DO 200 T=1 , 40 
FC ( I ) = FCU ) /MAX 
30n FWm=FW(l )/MAX 

600 Fnp.HAT{ 5H0C0C =G13.6, 1X6HC0C2 = G13.6 , 1X6HE0E = G13.6, 
11X7HTAUBAR= F5.1) 

C 

C TUT" T$ TO DtOT FC AND FW 

WRI T 5(inj t AUBAR,NTM,MAX t C0C2,E0E,FC,FW 
C 

C THERE CA RDS ARE FCR RESTART PROGRAM TR NEXT *** MARKS 

WRITE( 9) MOM t T,G,DVOL ,w,a ,vom,sf,co,to,ts,tf,tm ,dto,dtm,tn,max, 
1MAX6,C1,RH01, P, at ,ROT 1 , PDA, ABCC, R 
END FILE 9 
REWIND 9 

C *** END OF THE RESTART PROGRAM 

Rr T1JR M 
FMO 


SUBROUTINE X2D 1ST (TROT) 
COMMON /RANDOM/ R 
COMMON /IN 1X2/ ERO 
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c 

C in DETERMINE THE ROTATIONAL FREQUENCY(OM) from the 
C X2-D ISTR IBUTION ( K A T-SQUARE ) 

C 

C IROT= l»*SORT(,2‘5— 2. 48 503 EISNER DIALOG { 1-R ) )-.5 + .5 FOR 1 JUMP 

C I P. 0T= .50*SQRT(.25-2.4850_3El5*ER0*ALaG(l-R})-.25+.5 FOR EVEN JUMP 

C' I ROT = 2 #7 

C 10 R=RANF CO) 

C IROT= l.*SQRT C.25-2.48 503E15*ER0*AL0G( l.-R) ) -.1666 

C. RETURN 

C END 

*OECK SC A$$ 

FUNCTION SCATFRI VO, MOM) 

C IF (VO . GE iO . 0 ) SCATER=2.*AC0S{V0) 

C IF(VO.LT.O.O) SCATER=2.*{ 3. 1415926+ACOS (—VO ) ) 

SCATER=2.*AC0S(V0 ) 

RETURN 

END 

*DECK MOM$$ 

SUBROUTINE MOMENT 

COMMON /CV/ MAXtMAX6 t Cl»RH01 

COMMON /ANSWER/ MOM ,T ( 4) , G , DVOL , AT , PDA, R0T1 

COMMON /PART/ p{ 5001 ) 

C. 

C COMPUTES DESIRED MOMENTS 

C- 

DO POO 1=1 t MAX 
LU=5*T-3 
T ( 1 ) =T ( 1 J + l 
T(2)=T(?}+P(LU+3) 

T(3)=T(3)+P<LU+3)**2 
T ( 4 J = T f 4)+P { L U+4 ) **2+P ( LU+4J 
200 CONTINUE 
RE TURN 
END 


C 

C 

C SAMPLE INPUT DATA 
C 

If 

C F{ G/GM J-SELCT* ISELCT =1 , C=0 .4310, POO 




1000 


100 - 


. • 126 47? E- 2 43 632 *.03 3 .4651 9E- 22.44272E-141. 0 

3.0 p- 110.0 i. X1-. 

1-.XJ9- -F- 8-3;. 33 E-l 

»» 

ii 


46 



REFERENCES 


1. Deiwert, George S. : Reflection of a Shock Wave From a Thermally Accom- 

modating Wall; Molecular Simulation. Phys. Fluids, vol. 16, no. 8, 

Aug. 1973, pp. 1215-1219. 

2. Bird, G. A.: -Approach to Transitional Equilibrium in a Rigid Sphere Gas. 

Phys. Fluids, vol. 6, no. 10, Oct. 1963, pp. 1518-1519. 

3. Bird, G. A.: Direct Simulation Monte Carlo Method - Current Status and 

Methods, pp. 85-93; The Formation and Deflection of Shock Waves, 
pp. 301-311. In Rarefied Gas-Dynamics; Proc. Sixth Internatl. Symp., 
1968; Academic Press, N. Y., .1969. 

4. Deiwert, George S.; and Yoshikawa, K. K. ; Analysis of a Semiclassical 

Model for Rotational Transition Probabilities. Phys. Fluids, vol. 18, 
no. 9, Sept. 1975, pp. 1085-1088. 

5. Yoshikawa, K. K. ; and Itikawa, Y.: Monte Carlo Calculations of Diatomic 

Molecule Gas Flows including Rotational Mode Excitation. NASA TN D-8100, 
Jan. 1976. 

6. Bird, G. A.: Numerical Simulation and the Boltzmann Equation. In 

Rarefied Gas-Dynamics; Proc. Seventh Internatl. Symp., 1970. 

7. Hirschfelder, Joseph 0.; Curtiss, Charles F. ; and Bird, R. Byron: 

Molecular Theory of Gases and Liquids. John Wiley & Sons,. Inc., N. Y., 
1954. <- 

8. Parker, J. G. : Rotational and Vibrational Relaxation in Diatomic Gases. 

Phys. Fluids, vol. 2, no. 4, July-Aug. 1959, pp. 449-462. 

9. Lordi, John A.; and Mates, Robert E.: Rotational Relaxation in Nonpolar 

Diatomic Gases. Phys. Fluids, vol. 13, no. 2, Feb, 1970, pp. 291-308. 

10. Pearson, W. E. ; and Hansen, C. F. : Collision Induced Rotational 

Transition Probabilities in Diatomic Molecules. In Rarefied Gas- 
Dynamics; Proc. Eighth Internatl. Symp. 1972; Academic Press, N. Y., 

1974, pp. 167-175. 

11. Rabitz, H.: Rotational Energy Relaxation in Molecular Hydrogen, Jour. 

of Chem. Physics, vol. 63, no. 8, Oct. 1975. 

12. Yoshikawa, K. K. ; and Lee, K. L.: Collision Crossections for a Central 

Field Potential Function. ISAS RN-46, Univ. of Tokyo, 1977. 


47 



13. 


Vincenti, Walter G.; and Kruger, Charles H., Jr.: Introduction to 

Physical Gas Dynamics. John Wiley & Sons, Inc. , 1965. 

14. Rabitz, Herschel: Effective Potentials in Molecular Collisions. J. Chem. 

Phys., vol. 57-, no. 4, Aug. 15, 1972, pp. 1718-1725. 

15. Takayanagi, Kazuo: The Theory of Collisions Between Two Diatomic 

Molecules. ' Progr. Theor. Phys., vol. 11, no. 6, June 1954, pp. 557-594. 

16. Levine, R. D.; and Bernstein, R. B.: Molecular Reaction Dynamics. Oxford 

Univ. Press, 1974. 

17. Mark, Hans; and Paulissen, George T.: Electric Excitation of Certain 

Rare-Earth Nuclei by Protons. The Physical Review, vol. 100, no. 3, 
pp. 813-822, Nov. 1955. 

18. Polanyi, J. C.; and Woodall, K. B.: Mechanism of Rotational Relaxation. 

Jour, of Chem. Physics, vol. 56, no. 4, Feb. 1972. 


48 



TABLE I.- INITIAL DISTRIBUTIONS AND FLOW CORRELATION 


Figure 

Velocity* 2 

Rotational 0 " 

Energy- 

partition 

Temperature 

Remarks 

la & b 


f =f (T) 
r rm 

W 3/2 

T=T 

e 

Uniform flow 

2a & b 

W 

f ==frozen 
r 

— 

— 



3 

VW 

f H 
r rm 

E t /E r =3/2 

T =T 
t e 

Sound 

absorption 

experiment 

4 

f t- £ tm <T t ) 

f =f (T ) 
r rm r 

E t /E r »3/2 

T »T »T =0 
ter 

Shock wave 

5a & b 

f t- £ t» (T t> 

f =f (T ) 
r rm r 

E t /E r «3/2 

T. «T «T 
ter 

Free- jet 
expansion 

6a & b 

WV 

f H 

r rm 

E t /E r «3/2 

T. «T 
t e 

Chemical- 

fluorescence 

experiment 


a f tm(T) and frm(T) denote translational and rotational Maxwell-Boltzmann 
distributions with temperature T. 
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(a) Monte Carlo Results: ft(0,x) = Maxwellian where x = c/cg 

f r (0,j) = Boltzmann; T t = T r = 320 K. 


Figure 1.- Distribution functions for complete equilibrium. 
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(b) Time average of Monte Carlo results . 
Figure 1.- Continued. 
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(b) Time average of Monte Carlo results - Concluded. 
Figure 1,- Concluded. 
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(a) Delta function initial velocity distribution: ft-(0,x) = 1 

at x = /3/2. 


Figure 2 .- Monatomic gas simulation (rotational effect frozen). 
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(b) Double delta function initial velocity distributions; 
f t (0,xi) = 1/2 at xi = 1/2 
f t (0,x2) = 1/2 at X 2 = /TT/2. 


Figure 2 .- Continued 
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(b) Double delta function initial velocity distributions - Concluded. 

Figure 2.- Concluded. 
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Figure 3.- Maxwellian initial velocity distribution; delta function 
rotational energy distribution (f r (0,x) = 1 at j = 10) : 
Equipartition satisfied. 
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(a) Comparison with equilibrium distributions. 


Figure 5.- Maxwellian initial velocity distribution; Boltzmann 
rotational energy distribution (T t = 6 K, = 793 K) : 
Equipartition not satisfied. 
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(a) Comparison with equilibrium distributions - Concluded. 

Figure 5.- Continued. 
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(b) Comparison with local Maxwell-Boltzmann, distributions . 

Figure 5.- Continued. 
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Figure 5.- Continued. 
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(b) Comparison with local Maxwell-Boltzmann distributions - Concluded. 

Figure 5.- Concluded. 
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(b) Comparison with local Maxwell-Boltzmann distributions. 

Figure 6.- Continued. 
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(b) Comparison with local Maxwell-Boltzmann distributions 

Figure 6.- Concluded, 


Concluded. 
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